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I. INTRODUCTION 

Organic semiconducting materials have been the sub- 
ject of much attention in recent years due to the prospect 
of tailoring their photoelectric properties to specific 
applicationsP^ Doping organic semiconductors, by em- 
bedding electron- or hole-donor species in the material, 
provides a particularly effective way of tuning device 
properties.^ Dopant concentrations can be easily and pre- 
cisely controlled, although the effect on device perfor- 
mance is not always known a priori. 

Simulation is a powerful tool for understanding and 
optimizing device performance, allowing systematic stud- 
ies of device design that are experimentally infeasi- 
ble. Computational approaches taken by other au- 
thors include Monte Carlo simulation, 5 master equation 
approaches,^ and solution of the one-dimensional drift- 
diffusion-Poisson system of equations^ which is the ap- 
proach we follow here. This last approach is highly 
suited to the extraction of model parameters from ex- 
perimental data and device optimization, since it is less 
computationally-expensive than methods dealing with 
more dimensions. One-dimensional drift-diffusion mod- 
els originate from the theory of conventional crystalline 
semiconductors. However, when combining models for 
disordered materials with high doping intensities, numer- 
ical approaches originally developed for inorganic crys- 
talline materials may fail. We derive a generalization 
of existing methods, namely the Scharfetter-Gummel 
discretizatiorP and the Gummel-iteration map^, to over- 
come these problems and then compare our results with a 
state-of-the-art simulation approach presented in Ref . [TO] 
The numerical scheme that we develop here improves the 
numerical stability and computational efficiency of the 
approach. 



The rest of this paper is structured as follows: In 
Sec. |TlJ we introduce a model for steady-state charge 
transport in disorded semiconductors. In Sec. |III[ we de- 
scribe a scaling scheme to improve numerical behavior. 
We present in Sec. |IV| the proper Scharfetter-Gummel 
discretization of the model, accounting for the general- 
ized Einstein relation. In Sec. [Vj we derive the appropri- 
ately generalized Gummel iteration for the model. Nu- 
merical data for three example calculations are given in 
Sec. |VI| comparing the proposed method with that of 
Ref. 1101 Finally, we draw our conclusions in Sec. VII 
The procedure is written for hole-transporting and bipo- 
lar devices in Appendix [XJ and the scaling factors we use 
are tabulated in Appendix[B] Throughout this article, we 
use the symbol e for the elementary charge. 



II. MATHEMATICAL MODEL 

The starting point for the derivations of most macro- 
scopic charge transport models in semiconductors is the 
drift-diffusion-Poisson system introduced for inorganic 
crystalline semiconductors by van RoosbroeckJS] 
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Here, J is the particle current, p is the density, vb is the 
electrostatic potential, C is the dopant density (assumed 
to be positive in case of n-type doping) , fi is the electron 
mobility, and D is the diffusion constant. The subscripts 
n and p refer to electrons and holes, respectively, and the 
total charge current is given by — e( J n — J p ). To simplify 
the notation, we write the equations henceforth only for 
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electrons using the electron density n = p n . The follow- 
ing derivation is equally applicable for holes, since the 
equations in this case only differ by the sign of the drift 
term in the continuity equation and the carrier density 
in the Poisson equation. Appendix |A"| gives the extension 
of our results to bipolar and hole-transport devices. 

The model of Eqs. [I[j4] was derived for crystalline 
silicon-based semiconductors 1 ^! and has been applied in 
that area for many years with great success. Although 
there is no ab initio derivation, application of Eqs. [I]- 
[4] to organic semiconductors seems promising since the 
charge continuity and Poisson equations intuitively ap- 
ply in some form. Similar models are used for organic 
electronics in, e.g., Refs.HUlandHSl 

In fact, Eqs. [T}|4] must be extended in order to accu- 
rately describe disorded organic semiconductors.^ Avail- 
able models include the extended Gaussian disorder 
model (EGDM) and the extended correlated disorder 
model (ECDM)PEI] Both models assume a hopping 
mechanism for charge transport with a Gaussian density 
of states. The EGDM and ECDM introduce density- 
and field-dependence into the mobility and diffusion co- 
efficients /x(T,p,#) and D(T,/vH), respectively. The 



methods that we present in Sees. |III[ |IV| and |V| are appli- 
cable to arbitrary \x and D; here we discuss the EGDM 
and ECDM since they are recent and popular models. 

It should be noted that in the context of so-called deep 
traps, introduced e.g. in Refs. [10] and [TT1 the trapped 
carriers act as a stationary background charge and their 
treatment is completely analogous to the doping. Our 
method should therefore also lead to better performance 
in situations where deep traps are considered. 

We restrict the discussion to one-dimensional steady- 
state current simulations. It is useful to introduce the 
functions F\ and F%, writing the steady-state system as 
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where T is the temperature. When the system is ac- 
curately described by Maxwcll-Boltzmann statistics, the 
classical Einstein relation holds, 



D = V thf i , 



(7) 



with the thermal voltage Vth = fcsT/e. Note that in this 
equation the Boltzmann constant ks is taken to have 
units JK _1 . In order to align notation with other authors 
in the field, we will frequently use kg in units of eVK -1 . 
To still be able distinguish these two conventions we will 
use k in this case. 

When facing high charge-carrier densities, the effects 
of the Pauli exclusion principle become important and 
one is forced to resort to Fermi-Dirac statistics. In the 



context of organic semiconductors, it was first pointed 
out by Roichmann and Tessler^ that one should then 
use a generalized version of Eq. [7J 
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where tp n is the electron quasi-Fermi level (sometimes re- 
ferred to as the chemical potential; we use the term quasi- 
Fermi level to avoid confusion and to stress that it is a 
non-equilibrium quantity). Within the EGDM/ECDM 
framework, Eq. [8] is usually written as 



D = g 3 V th n , 
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and <?3 is termed the diffusion enhancement. We adopt 
this convention for clarity and compatibility with other 
authors. 

For a given density, the quasi-Fermi level is defined 
implicitly by the identity 
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where DOS(_E) represents the density of states. 

Equation [8] is called the generalized Einstein relation^ 
and accounts for the fact that charge carriers near the 
quasi-Fermi level are more likely to contribute to diffu- 
sion. The derivative occurring in Eq. [8] can be computed 
analytically by differentiating Eq. [10] The differentiated 
expression can be found in e.g., the Appendix of Ref. l20l 



III. SCALING 



In SI units, the absolute values of the carrier density 
[m~ 3 ] and electrostatic potential [V] can differ by over 
20 orders of magnitude. To avoid numerical problems we 
thus rescale the equations. Several scaling schemes have 
been proposed by other au thors for both inorganic crys- 
talline and organic devices J2U22I \y e £ n( j i^g^ the scaling 
leading to the most stable behavior is state-dependent. 

We suggest scaling the density by its maximum value 
(taking doping into account). We refer to the density 
scalingfactor as n sca i. Since we use a thermionic injection 
modeP^l, which can be coupled with barrier lowering^!, 
the boundary values and therefore maxima of the density 
are not known a priori. The scaling is therefore applied 
at every iteration. The potential is scaled by V'scaij the 
maximum value of (the absolute value of) the applied 
voltage V app i and the thermal voltage. This ensures that 
the scaled potential is reasonably small even for high ap- 
plied voltages. On the other hand, we force the scaling 
parameter to be at least the thermal voltage to avoid 
scaling by zero. Since the mobility may differ strongly 
between materials, we scale all mobilities by their zero- 



3 



field, zero-density limit fig. Finally, space is scaled by 
L, the total thickness of the device. The full list of the 
scaling parameters including the resulting scaling of the 
current density is given in Table [n] 

The scaled system of equations reads (using the same 
symbols for the scaled quantities) 
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Here, A is the scaled Debye length. 

This scaling ensures that the potential and density vary 
between and 1 for most physical situations. Failing to 
implement a suitable scaling scheme can result in nu- 
merical instability and loss of accuracy. Furthermore, it 
reduces the number of constants used in the actual com- 
putations. Except where explicitly stated, we henceforth 
use the scaled quantities. 



IV. DISCRETIZATION 

To solve Eqs. [5] and [6] (note that these are the unsealed 
equations), Scharfetter and GummeP proposed a scheme 
for the case of constant mobility and diffusion. Instead 
of using central or upwind differences to approximate the 
current, a weighted difference depending on the field is 
used. It can be derived as a solution of the boundary 
value problem (BVP) 



J i+ i = A* (nVip - VthVn) 
n(x = Xi) = Hi , 
n(x = Xi+i) = n i+1 , 
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which describes the physics on a sub-interval for given 
carrier densities at the boundaries. The BVP of Eq. [14] 
is first order, but we can prescribe two boundary val- 
ues because Jj+i is a free parameter of the system, and 
there is only one value of J i+ i that admits a solution. 
This value can be used as an approximation of the cur- 
rent, assuming that the mobility, diffusion coefficient and 
electric field are locally constant. Expressed in terms of 
11 and ip, it reads 
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The scheme can be viewed as an upwind scheme, where 
the amount of upwind difference depends on the local 
strength of the field. This is easily seen by considering 
the asymptotic behavior of when the field is large 
or small, 
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so that for small fields we use a classic central difference 
approximation for the gradient of the density, i.e., the 
diffusive part of the current. When the magnitude of the 
field is large, we obtain 
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which are exactly upwind approximations for the density. 

Upwind schemes introduce artificial diffusion into the 
solution) 2 ^ which is a source of error. In the case of con- 
stant mobility and the simple Einstein relation of Eq.[7 it 
can be proven that the discretization introduced in Eq. 15] 
yields the optimal artificial diffusion (i.e., the minimum 
required for numerical stability).^ 

Within high-density capable models such as 
EGDM/ECDM, the generalized Einstein relation of 
Eq. [8] increases physical diffusion. Use of Eq. [15] to 
approximate the current is then suboptimal in terms 
of the accuracy of the solution because more artificial 
diffusion is added than is required for stability. Using 
the generalized Einstein relation Eq. [S] and applying the 
scaling introduced in Sec. |III| we obtain the following 
generalization of the BVP of Eq. PHI 
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This corresponds to the current in Eq. [TT] 

Assuming n, g 3 and V^> to be constant on each interval, 
or put equivalentl y, ap proximating them by their value at 



1, we solve Eq. 19 as before to get a properly-adjusted 



discretization scheme 
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This scheme accounts for the additional physical diffu- 
sion in our system by decreasing the artificial diffusion, 
improving the accuracy of the current. 
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Figure 1. The arrows /* and 4- represent solution of the 
Poisson equation and solution of the steady-state continuity 
equation, respectively. Dashed arrows show a quantity being 
carried forward for use in the next iteration. 



In the EGDM/ECDM, \i = n {T)g 1 (n,T)g 2 (Vip,T). 
In this case gi and 53 at half-integer gridpoints can be 
obtained by various averaging methods. Throughout this 
work we use the geometric average 
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which is an exact interpolant for exponentially varying 
functions. Note that since 172 depends on the gradient 
of the potential, it is already known on half-integer grid- 
points only. 



V. ITERATIVE SOLUTION OF THE DISCRETIZED 
SYSTEM 

To solve the system of discretized equations, Eqs. [5] 
and[6j using the discretization of Sec. |IV| one can either 
use the Newton algorithm or the so-called Gummel iter- 
ation, a physically-motivated variant of the Gauss-Seidel 
relaxation schemeP 

Numerical comparisons as in the work of Knapp et 
alW^ and additional theoretical considerations show that 
Newton's method converges in fewer iterations than the 
Gummel method. However, Newton iterations are gener- 
ally more expensive to compute. Furthermore, Newton's 
method often requires initial values in the vicinity of the 
solution in order to obtain convergence. 

This is a severe problem for optimization and parame- 
ter extraction, where the solver is routinely called several 
hundred times. In these cases a single failure to con- 
verge can be catastrophic for the optimizer. We there- 
fore derive an iterative solution scheme in the spirit of 
GummelP As a starting point we will use the continuous 
(non-discretized) problem to avoid working in spaces of 
grid functions. 

The idea is to alternately solve F\ for n at a given 
ip and F2 for ip and given n until self-consistency is 
achieved. 

Since the drift-diffusion-Poisson system of equations 
can be nonlinear (even when /i and D are constant; there 
is implicit dependency of e.g., ip on n), we must linearize 
first. Our iteration starts with the Poisson equation. 



A. Linearization of the Poisson equation 



The essence of the Gummel iteration is to consider the 
implicit nonlinearity of the Poisson equation explicitly, 
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In the original derivation, the nonlinearity of n in ip is 
revealed by the explicit use of Maxwell-Boltzmann statis- 
tics. This does not apply in our case, and we must use 
Eq. 10 instead. However, since even in Eq. 10 the de- 



pendence on ip is not explicit, we must first introduce 
another quantity, the quasielectrochemical potential £ n , 
as defined in e.g., Ref. 



Cn ■= ~ # 
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The quantity Cn can be interpreted as the single driv- 
ing potential for the charge carrier transport. One must 
be aware that there are several conventions in defining 
these levels, which are not consistent with each other. 
We use the distinction between the different levels as in- 
troduced in Ref. However, we use a definition of the 
quasi-Fermi level consistent with Ref. [201 since the 53- 
prefactor that we use in our computations is introduced 
there. Therefore, the quasi-Fermi level and quasielectro- 
chemical potential are switched with respect to those of 
Ref. [251 The difference stems from the choice of refer- 
ence point with respect to which we define the DOS. We 
need here to reformulate the Gauss-Fermi integral of 
Eq. [lOj Since this equation uses an unsealed ip we need 
to unscale our variable when we put it in 
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where A^ites is the site density. We deliberately use q here 
instead of e because its numerical value is 1. This can 
be understood from the consideration that all energies in 
the numerator of the exponential in Eq. [24] are in units of 
eV, whereas the product eip sca \ip has unit J. These units 
differ by the numerical value of the elementary charge e, 
hence q in Eq. [23] is 1C and can be omitted in numerical 
calculations. 

We find it useful to define the 'Poisson equation oper- 
ator' 



Poiss (ip, n(ip)) := AiP - A" 2 (n(^) - C) , 



(25) 



which we wish to linearize at the current iterate 
(tp^, n^), where k denotes iteration number. To avoid 
unnecessarily lengthy formulas, we introduce the nota- 
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The linearization is performed using functional 
derivatives,^ 
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cation operator. Since the central difference approxima- 
tion of the Laplacian and the pointwise approximation of 

are linear and continuous as projections onto finite 

dimensional spaces, the discretization commutes with the 
linearization. 
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g 3 factor. Comparing Eqs. [T0|and[24} we get 
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Equation 30 is highly convenient, since 53 and n are com- 
puted elsewhere and can be reused here at virtually no 
cost. The linearization process thus doesn't require com- 
putation of any additional quantities. 



Using Eq. 30 we get for the derivative of the Poisson- 
operator 
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To com plete the linearization, we insert this result into 



Eq. 27 set Poiss('0, n{ip)) = and -0 = ip^ k+1 \ and solve 



for This yields 
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We now switch back to the discretized system. As men- 
tioned before, the order in which we linearize and dis- 
cretize does not matter. We discretize using the following 
expressions, 
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and solve for ip( k+1 \ setting the first and the last (spa- 
tial grid) values of the quantity in the square brackets of 
Eq. [32] to the boundary values. The boundary values of 
ip and n are determined by the applied voltage and the 
charge injection, respectively. The matrix elements not 
shown in Eq. [34] are zero. 

Equation [32] which we derived by linearizing the pois- 
son operator, properly generalizes the Gummel iteration 
scheme to the case of a generalized Einstein relation. 



B. Linearized Continuity equation 

Since it is comparatively straightforward, we give the 
linearization for the continuity equation on the dis- 
cretized level. We do not consider the implicit nonlinear- 
ity via ip( n ) here, and linearize the mobility by assuming 
that ^ k+1 \ the mobility at the fc+l-th iteration, depends 
on and not n' fe+1 - 1 . Before giving the full expression 
for the linearized and discretized current, it is convenient 
to introduce the function 
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Substituting this definition into Eq. 20 the discretized 
and linearized current is given by 
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Figure 2. (Color online) Simulation results with the proposed algorithm. Where the method of Knapp et al. (not shown) 
converges, the results coincide with ours. The applied voltage in the simulation was set to 3V. Model parameters are given in 
Table □ 



This expression is linear in n^ k+1 K Therefore 
the central approximation of the continuity equation 



0, which is what is solve d in p ractice, is 



define the 
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clearly also linear in n^ k+1 \ Equations |32f[36 
discretized and linearized system and are accessible to 
direct solution methods. 

A diagram illustrating the iteration scheme resulting 
from the linearizations presented in the preceding sec- 
tions can be found in Fig. [T] 



VI. SIMULATION RESULTS 

To highlight the advantages of the schemes described 
above, we have simulated a simple single-layer structure 
using the ECDM model. We consider three different situ- 
ations: First, we simulate a low density situation without 
doping — here we show that our approach reproduces 
known results with a small gain in computational effi- 
ciency. Second, we simulate a device with high carrier 
densities due to a low injection barrier — here we ob- 
serve a significant performance gain with the proposed 
method. Finally, we simulate a strongly-doped device, 
with approximately every 100th molecule replaced by a 
dopant. In this third case, conventional methods fail but 
our procedure converges. The parameters used for the 
simulations are collected in Table |T] The resulting charge 
carrier densities and potentials at 3V are shown in Fig. [5] 
All of the device simulations, including those using the 
conventional iteration, used the discretization described 
in Sec. IIVI and had identical initial values. Our calcula- 



Table I. Simulation parameters for the three example cases 
that we consider. Simulation results are given in Figs. [2}j4] 



Parameter 


(1) 


(2) 


(3) 


Electrode WF [eV] 
Doping intensity 


-5 
0% 


-4.6 
0% 


-4.6 
1% 


Mobility model 
Injection model 
LUMO [eV] 
Site density [m -3 ] 
DOS width [eV] 
Temperature [K] 
Device Length [nm] 
ECDM-C parameter 
/i [m^Vs)- 1 ] 


ECDM 
Thermionic 
-4.5 
2 x 10 27 
0.13 
300 
100 
0.29 
4.5 x 10~ 6 



tions did not include the barrier-lowering effect of image- 
charges in the metal contacts. Its inclusion would have 
the effect of increasing the carrier injection, further en- 
hancing the computational advantage of our method. 



A. Convergence 

In Fig. [3j we compare the behavior of the method pro- 
posed in this paper, with the method presented by Knapp 
et al, 10 which we call the 'conventional method'. 

The main difference between the approaches lies in the 
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Figure 3. (Color online) Plot of convergence speed in Z 2 -stepsize of the carrier density for the three different scenarios given 
in Table [I] The 'proposed' method is that described in this paper, while the 'conventional' method is the iteration given by 



Knapp et aZp^ Parameter values are given in Table [I] 



linearization of the Poisson equation. Knapp et al. use 
the expression (written in the framework of our scaling 
for the reader's convenience) 
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instead of Eq. 32 This differs from Eq. 32 by the absence 
of <?3 in the denominators and the doping profile, which 
wasn't considered in Ref. [TOl The effect can be viewed 
as a different damping. However, one should keep in 
mind that our scheme is not derived from a perspective 
of damping but from the use of functional derivatives for 
the Poisson operator introduced in Eq. [25] 

A similar difference can be seen in the discretization of 
the continuity equation, where an additional (73 turns up 
in the denominator of the argument of the exponential 
function. As mentioned above, this prevents the intro- 
duction of unnecessary artificial diffusion. 

In all three device simulations, the proposed method 
is superior to the conventional method in terms of con- 
vergence speed, with the difference being greater when 
the carrier densities are high. Our convergence criterion 
was that the / 2 -stepsize in the scaled charge density fall 
below 1CT 7 . 

We find that the stability of the generalized iteration in 
high-concentration regimes is greatly enhanced in com- 
parison with the conventional iteration. For example, the 



conventional approach fails to converge for doping inten- 
sities > 1%, whereas our method converges for arbitrary 
doping intensities. The benefits of properly incorporat- 
ing diffusion into the numerical scheme could also be im- 
portant for the simulation of multilayer devices, and in 
particular those with doped layers. 

Even in the absence of doping, for large values of charge 
carrier injection a speed-up of convergence is observed. 
In the low density limit, the method blends into the one 
presented in Ref. fTTJl because g$ — > 1, as n — > 0. 



B. Grid convergence 

Since we introduced a scheme using upwind stabiliza- 
tion, it is interesting to study the effect of varying grid- 
point density on the behavior of the solution. We have 
studied the effect of the number of gridpoints on the cur- 
rent density for the three cases shown in Table [TJ We 
investigate uniform grids with up to 1000 gridpoints. As 
shown in Fig. [4j we observe convergence using the pro- 
posed method in all three cases. However, for high car- 
rier densities the relative error in the current oscillates 
slightly with varying grid resolution. This is likely to be 
a result of the averaging procedure used on the mobility 
prefactors 51-3, and is a topic for further investigation. 
The oscillations are, however, smaller than the general 
trend of convergence. 
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Figure 4. Convergence of the calculated current as the numer- 
ical grid is refined. Parameter values are given in Table [I] The 
relative error is defined as <5 J rG i = ( Jjv — J 10 3 ) / J 10 3 , where Jjv 
is the calculated current with N grid points. 



VII. SUMMARY AND CONCLUSIONS 



We have introduced a numerical method consisting of a 
scaling scheme, a discretization scheme, and a fixed point 
iteration for the van Roosbroeck drift-diffusion-Poisson 
system with arbitrary density- and field-dependent mo- 
bility functions. The method properly accounts for the 
implications of Fermi-Dirac statistics by incorporating 
the generalized Einstein relation. The iteration is de- 
rived by linearizing the Poisson equation with explicit 
treatment of the nonlinear dependence of the density on 
the potential. 

Failure to take consequences of the generalized Ein- 
stein relation into account in the iteration results in nu- 
merical instability when the density is large, e.g., in the 
case of strong doping or large trap density, whereas the 
generalized method converges for arbitrary doping, fn 
other cases, we find that the generalized method is sig- 
nificantly more efficient than the conventional method 
while reproducing converged current values. 



Appendix A: Numerical Scheme for positive carriers 

While the modelling details of bipolar devices (charge 
generation/recombination) are outside the scope of this 
paper, we show how to extend the method presented 
above to hole-transporting or bipolar devices. For the it- 
erative procedure, we include a hole continuity equation, 
which has the same mathematical form as the electron 
continuity equation, except the sign of the drift term is 
changed. As for electrons, we use Scharfetter-Gummel 
for the discretization and linearize the mobility by using 
the values of the last iteration. The analogue of Eq. [36] 
for holes is 



j(fe+i) 



' (fc+1) n {k) (fc+1)' 
Pi ~ U p,iPi+l 



' , (fc+1) ,(fc+l)' 



I - G 



(fc) 



(Af) 



This allows us to solve the for the new iterate p( k+l ^ in 
the same way as we did for the electron density in Sec.|V] 
In the linearization of the Poisson equation, we do not 
only need to linearize n(ip) but in an analogous manner 
also p(ip). The effect of the potential on the two carrier 
densities is the same except for a reversed sign, i.e., the 
hole analogue of Eq. [28] has a minus sign. However, since 
n and p appear in the poisson equation with different 
signs, they appear in the linearized iteration in the pref- 
actor multiplying ip^ on the same footing (i.e., + x — 
or — x +). We obtain 



(fc+1) 



d Poiss 



dtp 



A 2 V 



,(fc) _ p (fc) 



C 



p r n 



(fc) 



A 2 « 

i/3,n 



p r p^ 



(fc) 



(A2) 



Here, the linearized Poisson operator reads 



d Poiss 



dip 



p r n 



(fc) 



A 2 « 

y3,n 



,(fc) 

~Jkj 
93,P 



Pr P y 

A 2 J fc ) 



(A3) 



The discretization and solution of the Poisson equation 
then works exactly as described in Sec |V| 
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Appendix B: Scaling factors 

Here we list the model scaling factors, which we use 
to force densities and potential to vary between zero and 
one. Note that the density and current scaling factors 
have the potential to change during the calculation, so 
should be updated after every iteration. 
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Table II. Scaling factors for quantities considered in the 
model. The symbols in the second column are defined as 
the quantities in the third column. 



12 A. Juengel 
Chap. 5, pp 



Variable Symbol 



Scaling factor 



Potential 
Density 

Length 

Mobility 

Current 



V'scal 
^scal 

L 

"seal 



max{|Va PP i|, V t h} 
max{p(x = 0), p(x = i), 
max^^L] C(x)} 
L 
Mo(T) 

/^scal^scal^scal-^ 



REFERENCES 

1 A. Pron, P. Gawrys, M. Zagorska, D. Djurado, and R. De- 

madrille, |Chem. Soc. Rev. 39, 2577 (2010)| 
2 S. van Mensfoort, J. Billen, M. Carvelli, S. Vulto, R. Jansscn, 

and R. Coehoorn, J. Appl. Phys. 109, 064502 (2011). 
3 J. E. An thony, A. F acchetti, M. Heeney, S. R. Marder, and 

X. Zhan, [Adv~ Mater. 22, 3876 (2010)) 
4 C. Chan, W. Zhao, A. Kahn, and I. Hill, Appl. Phys. Lett. 94, 

203306 (2009). 

5 J. van der Hoist, F. van Oost, R. Coehoorn, and P. Bobbert, 
Phys. Rev. B 83, 085206 (2011). 

6 S. van Mensfoort, V. Shabro, R. de Vries, R. Janssen, and R. Co- 
ehoorn, J. Appl. Phys. 107, 113710 (2010). 

7 J. Cottaar, R. Coehoorn, and P. Bobbert, Org. Elec. 13, 667 
(2012). 

8 D. Scharfetter and H. Gummel, IEEE T. Electron. Dev. 16, 64 
(1969). 

9 H. Gummel, IEEE T. Electron. Dev. 1, 455 (1964). 
10 E. Knapp, R. Hausermann, H. Schwarzenbach, and B. Ruh- 

staller, J. Appl. Phys. 108, 054504 (2010). 
n W. van Roosbroeck, Bell System Tech.J, 29, 560 (1950). 



"Drift-diffusion equations," (Springer, 2009) 
9-127. 

13 Y. Xu, T. Minari, K. Tsukagoshi, R. Gwoziecki, R. Coppard, 
M. Benwadih, J. Chroboczek, F. Balestra, and G. Ghibaudo, 



|Journal of App lied Physi cs 110, 01 4510 (2011) 

1<*W «^V, OT 1\/T T? Qkl^i^ov onH T T Ron^W 



Phys. Today 44,| 



4 H. Scher, M. F. Shlesinger, and J. T. Bendler, 
| 26 (1991)1 

1J R. Coe hoorn, W. F. Pasveer, P. A. Bo bbert, and M. A. J. 
Michels, |Phys. Rev. B 72, 115206 (2005)| 

16 M. Bouhassoune, S. van Mensfoort, P. Bobbert, and R. Co- 
ehoorn, Org. Elec. 10, 437 (2009). 

17 M. Schober, M. Anderson, M. Thomschke, J. Widmer, M. Furno, 
R. Scholz, B. Liissem, and K. Leo, Phys. Rev. B 84, 165326 
(2011). 

18 Y. Roichman and N. Tessler, Appl. Phys. Lett. 80, 1948 (2002). 
19 N. W. Ashcroft and D. N. Mermin, Solid State Physics (Saunders 

College Publishing, 1976). 
20 S. L. M. van Mensfoort and R. Coehoorn, |Phys. Rev. B 78, | 
| 085207 (2008)| 

•"F. Brezzi, L. Marini, S. Michelettio, P. Pietra, R. Sacco, and 
S. Wang, "Handbook of numerical analysis, vol. xiii (numerical 
methods for electrodynamic problems)," (Elsevier Science, 2005) 
Chap. Finite element and finite volume discretizations of Drift- 
Diffusion type fluid models for semiconductors. 

22 J. Bonham and D. Jarvis, A ust. J. Chem. 30, 705 (1977). 

23 J. Scot t and G. G. Malliaras, |Chemical Physics Letters 299, 115| 

| (1999)| 

^P. R. Emt age and J. J. Q'Dwyer, | Phys. Rev. Lett. 16, 356 (1966) 

25 G. Smith, Numerical Solution of Partial Differential Equations: 
Finite Difference Methods . Oxford Applied Mathematics And 
Computing Science Series (Clarendon Press, 1985). 

26 E. Doolan, J. Miller, and W. Schilders, Uniform numerical meth- 
ods for problems with initial and boundary layers Advances in 
numerical computation series (Boole Press, 1980). 

27 E. Kna pp and B. Ruhstaller, |Opt. Quant. Electron 42, 667] 

| (2011)| 

zo C. F. W oellner and J. A. Freire, |J. Chem. Phys. 134, 084112| 
| (2011)| 

2a R. Courant and D. Hilbert, "The calculus of variations," in Meth- 
| ods of Mathematical Physics\ (Wiley- VCH Verlag GmbH, 2007) 



pp. 164-274. 



